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Abstract: We compute the electromagnetic form factor of a "pion" with mass m n = 

330 MeV at low values of Q 2 = —q 2 , where q is the momentum transfer. The computations 

^^ ' are performed in a lattice simulation using an ensemble of the RBC/UKQCD collabora- 

tion's gauge configurations with Domain Wall Fermions and the Iwasaki gauge action with 



an inverse lattice spacing of 1.73(3) GeV. In order to be able to reach low momentum 
transfers we use partially twisted boundary conditions using the techniques we have de- 
veloped and tested earlier. For the pion of mass 330 MeV we find a charge radius given 
by (7^)330 MeV = 0.354(31) fm 2 which, using NLO SU(2) chiral perturbation theory, ex- 
trapolates to a value of (r 2 ) = 0.418(31) fm for a physical pion, in agreement with the 
experimentally determined result. We confirm that there is a significant reduction in com- 
putational cost when using propagators computed from a single time-slice stochastic source 
compared to using those with a point source; for m^ = 330 MeV and volume (2.74 fm) 3 we 
find the reduction is approximately a factor of 12. 

Keywords: Lattice QCD, Nonperturbative Effects, Pion Physics, Electromagnetic 
Processes and Properties . 
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1. Introduction 

In this paper we compute the electromagnetic form factor of a "pion" with mass m n = 
330 MeV at low values of Q 2 = —q 2 , where q is the momentum transfer. The computations 
are performed in a lattice simulation using an ensemble of the RBC/UKQCD collabora- 
tion's gauge configurations with Domain Wall Fermions and the Iwasaki gauge action with 
an inverse lattice spacing of 1.73(3) GeV (see section 3.2 for brief details of the simulation 
and ref. [1] for a full discussion 1 ). The action has good chiral and flavour symmetries and 
as demonstrated in ref. [1] a mass of 330 MeV is well within the regime where NLO SU(2) 
chiral perturbation theory holds for other physical quantities, such as the meson masses, 
decay constants and the kaon's bag parameter. 

In order to be able to reach low momentum transfers we employ partially twisted 
boundary conditions using the techniques developed and tested in ref. [2]. Previous lat- 
tice computations have used quarks satisfying periodic boundary conditions and therefore 
obtained form factors at much larger values of Q 2 (see however, the preliminary study 



1 In ref. [1] the pion mass corresponding to the bare quark mass used in the present study was found to 
be 331(6) MeV, where the error is dominated by the uncertainty in the lattice spacing. In the text we refer 
to this meson as having a mass of 330 MeV, while in the analysis we treat the fluctuations in the mass using 
a jackknife procedure. 



with twisted boundary conditions presented by the European Twisted Mass Collabora- 
tion (ETMC) [3]). For the pion with m n = 330 MeV we find for the charge radius, 
(?^)330 MeV = 0.354(31) fm 2 . We then use NLO chiral perturbation theory to obtain the 
form factor and charge radius of a physical pion, finding 

(r 2 ) =0.418(31) fm 2 , (1.1) 

in agreement with the experimentally determined value. 

The power of the technique is demonstrated in fig. 2 where the data points are obtained 
from our simulation. The dashed vertical line is the minimum value of Q 2 (Q^m) which is 
accessible with periodic boundary conditions. From the figure we see that the form factor 
can be obtained at arbitrarily small values of Q 2 and also that the results obtained with 
twisted boundary conditions join smoothly onto those obtained by performing the Fourier 
sum in the conventional way (i.e. onto the data point on the dashed line). In this paper we 
focus on the pion's electromagnetic form factor, but we anticipate that the technique used 
here will also have important applications to the calculation of other flavour non-singlet 
form factors at arbitrary values of momentum transfer, such as those which appear in K^ 
semileptonic decays [2]. 

For this calculation we use propagators generated from a single time-slice stochastic 
source in addition to standard point source propagators. We compare the cost, obtaining 
similar errors for the pion mass, the normalization constant of the vector current, Zy, and 
the pion's electromagnetic form factor at Q 2 ^^, finding, for m^ = 330 MeV and volume 
(2.74 fm) 3 , a gain of approximately a factor of 12 in favour of the noise source propagators. 
A gain was also found in the preliminary study presented by the ETMC collaboration in [3] 
and in the recent publication by the UKQCD collaboration [4]. 

The plan for the remainder of this paper is as follows. In the next section we briefly 
review the use of partially twisted boundary conditions to compute the form factor at 
values of Q 2 which are inaccessible with periodic boundary conditions [2]. The details 
of the computation, the parameters of the simulation and the results for the form factor 
for the 330 MeV pion are presented in section 3. The use of noise source propagators to 
evaluate the three point correlation functions from which the form factor is obtained is 
briefly described in section 3.3 and a comparison of the relative cost of using point source 
and noise source propagators to obtain results with the same statistical error is given in 
section 3.5. The use of NLO chiral perturbation theory to obtain the form factor and 
charge radius for a physical pion from that with mass 330 MeV is described and performed 
in section 4. Finally, in section 5 we present our conclusions. 



2. Twisted boundary conditions and the form factor at small Q 2 

The electromagnetic form factor of the pion, f' K1T {q 2 )-, is defined by the matrix element 

(vt+CpOI^K+Cp)) = (p+J>%rV), wher e q 2 = -Q 2 = (p - p'f (2.1) 
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Figure 1: Sketch of the valence quark flow in the electromagnetic form factor of the pion. There 
is a similar contribution in which the current is on the antiquark line and the spectator is a quark. 

and V^ = ■^u'y^u — -^dj^d is the electromagnetic current. In a finite volume with periodic 
boundary conditions for the quark fields, the accessible pion momenta are given by 

p=(E n ,p H ) = {E n ,(2ir/L)n) and p> = (E n ,,p H r) = (E n r, (2n/L) ft') (2.2) 

where n and n' are vectors of integers, L is the spatial extent of the lattice and E n and E n i 
are the corresponding energies (E 2 = m 2 + (2tt/L) 2 \n\ 2 and E 2 , = m 2 + (2ir/L) 2 \n' | 2 , 
where m^ is the mass of the pion), so that q 2 can only take the corresponding discrete values. 
In particular the minimum non-zero value of Q 2 is given by Q^m = 2m 7T (y / m 2 + (2tt/L) 2 — 
m n ), which for the parameters of our simulation is about 0.15 GeV 2 . In this paper we 
study the form factor at small Q 2 (and in particular for Q 2 <C 0.15 GeV 2 ), using the new 
technique proposed in [2] which allows one to carry out lattice computations at arbitrarily 
small values of Q 2 . We now briefly review this technique. 

In order to reach small momentum transfers, we use partially twisted boundary condi- 
tions [5,6], combining gauge field configurations generated with sea quarks obeying periodic 
boundary conditions with valence quarks with twisted boundary conditions [5-14]. The 
valence quarks satisfy 

q(x k + L) = e ie «q(x k ), (A; = 1,2,3), (2.3) 

where q represents one of the degenerate up or down quarks. We have demonstrated in 
section 2.3 of ref. [2] that it is possible to introduce twisted boundary conditions indepen- 
dently for the three valence quarks and antiquarks, i.e. 9\ for qx, 62 for q^ and #3 for qj, 
in fig. 1. In our study it will be sufficient to set #3 = so that the spectator quark or 
antiquark satisfies periodic boundary conditions. By varying 9\ and 62 we are able to tune 
the momenta of the initial and final pions continuously. 

The dispersion relation for a meson with twisting angle 6 takes the form [8, 10], 



E^ = ^ml+[p H + -\ , (2.4) 

where m v is the pion mass and p^ is the meson momentum induced by Fourier summation. 
For the matrix element in (2.1) with the initial and the final meson carrying momenta 
P '= PH + 0/L and p' = pfi> + O'/L respectively (where 6 = 9\ — 63 and 9' = 62 — 0$), the 
momentum transfer between the initial and the final state meson is 

q 2 = (p- p'f = (E^p) - E^p' )) 2 - ((p H + 9/L) - (p H , + O'/L)) 2 . (2.5) 



The finite-volume corrections with partially twisted boundary conditions decrease ex- 
ponentially with L similarly to those with periodic boundary conditions [5] . 

3. The computation and results 

In this section we present the details of the computation of the electromagnetic form factor 
of a pion with mass m^ = 330 MeV. In the first subsection we explain which correlation 
functions are computed in order to be able to extract the form factor. The parameters 
of the simulations are presented in subsection 3.2 and a brief introduction to the use of 
noise-source propagators is given in subsection 3.3. In section 3.4 we present our results 
for the form factor. Finally in section 3.5 we compare the computational cost of computing 
correlation functions with point source and Z(2)-wall source propagators. 

3.1 Correlation functions 

In order to determine the form factors we compute two- and three-point correlation func- 
tions. The two-point function is defined by 

C n (t,p) = Y / e^ £ (O n (t,x)Oi(0,6)) = JJ^- ( e -^(p)* + e -^(p)(r-')) , ( 3 .1) 

X 

where O n = dj^u is a local pseudoscalar interpolating operator for the pion. We have 
assumed that t and T — t (where T is the temporal extent of the lattice) are sufficiently 
large for the correlation function to be dominated by the lightest state (i.e. the pion). The 
constant Z w is given by Z^ = ( it | Oj(0, 0) | ) . The three-point function is defined by 

C^(t,t f ,p,p') = Z V ^j?<*f-QjP*(O v (t f ,a! f )V A (t,x)OU0,d)) 

X},X 

-4^Sk < ^' )|Vi(0,l ^ ,) 

x \e{tf - t) e -EAp)t^(p')(t f -t) _ 0{t _ tf) e -E^(p)(T-t)-E n (p')(t-t f )j ^ (32) 

where V4 is the time component of the bare electromagnetic current and where, without 
loss of generality, we have placed the source at the origin. Again we assume that all the 
time intervals in (3.2) are sufficiently large for the lightest hadrons to give the dominant 
contribution. As explained in the following paragraph, Zy is the normalization factor by 
which the bare lattice current needs to be multiplied in order to obtain the physical current. 
The normalization factor Zy can readily be obtained as follows. For illustration we 
take < t < tj < T/2, in which case Zy is defined by 

Zy= p^L . (3.3) 

eg, («,*/, 0,0) 

In the numerator we use the function C n (t,p) = C\(t,p)— 2E <v) e _jE,r ^( T ~^ where Z n and 
E-wip) = \/fn^ + p 2 are determined from fits to C w (t, ). For tf < T/2 this proves to be an 



effective and numerically stable procedure for the subtraction of the contribution from the 
backward propagating meson to C w (tf,0) in the numerator of (3.3). (For tf = T/2 it is 
natural instead to use CV(t,p) = ^CV(t,p) in (3.3).) The superscript B in the denominator 
indicates that we take the bare (unrenormalized) current in the three-point function. 

In the following subsection we introduce the three datasets which we use for our anal- 
ysis. For data set A we do not use twisted boundary conditions, setting p' = and 
determining the pion form factor from the ratio of correlation functions 

C n7T {t,t f , 0,0) C w (t,p) 
For data sets B and C we use 



IC nn (t,tf,p,p')C n7T (t,t f ,p',p) 



2VE n (p)E n (p') X r y\ + -XT in — r{d 2 ){E«(p) + EAp')) , (3.5) 

V C 7T {t f ,p)C 7T {t f ,p') 

(called ratio R\ in [2]). Both ratios approach a constant for sufficiently large time intervals. 

3.2 Parameters of the simulation 

The computations described in this paper were performed using the ensemble with light 
quark mass am u = am^ = 0.005 and strange quark mass am s = 0.04 from the set of 
Nf = 2 + 1 flavour Domain Wall Fermion [15-17] configurations with (L/a) 3 x T/a x L s = 
24 3 x 64 x 16 which were jointly generated by the UKQCD/RBC collaborations using 
the QCDOC computer [18-21]. The bulk of the correlation functions were evaluated on 
the UK Research Councils' HECToR Cray XT4 computer, with the set completed using 
a University of Edinburgh BlueGene/L system as well as QCDOC. A detailed study of 
the light-hadron spectrum and other hadronic quantities using these configurations has 
recently been reported in ref. [1]. The gauge configurations were generated with the Iwasaki 
gauge action [22, 23] at (3 = 2.13 corresponding to an inverse lattice spacing of a -1 = 
1.729(28) GeV. The resulting pion mass is m n ~ 330 MeV. In our numerical evaluations 
we use the masses measured directly on our data sets, which can be determined from the 
entries in tab. 2 and are fully consistent with the value reported in [1]. We use the jackknife 
technique to estimate the statistical errors. 

In the following we distinguish three sets of correlation functions as specified in tab. 1. 
Set A was generated with point sources and sinks. We started the measurements for three 
different source positions on trajectories 900, 905 and 910, respectively, measuring on every 
40th trajectory in each case and averaging three consecutive measurements over the sources 
into one bin. The initial pion carries momentum \p\ = 0, -£- or v2-f and the final pion is 
at rest. For this dataset we do not use twisted boundary conditions at all. 

For data sets B and C we used a Z(2) x Z(2) noise wall source as explained in sec- 
tion 3.3 and a point sink. For data set B we started the measurement chains for the eight 
source positions specified in tab. 1 on trajectories 900, 905, 910, . . . , 935. Data set C is a 
subset of set B which starts with four source positions on trajectories 900, 910, 920 and 930, 
respectively. In each case we measured on every 40th trajectory and averaged the correla- 
tion functions over the chains into bins of eight and four measurements, respectively. The 



set trajectories on £ src = A N n 



A 


900-4460 


20 


537 


0,16,32 


B 


1000-6840 


40 


1176 


0, 54, 20, 14, 56, 26, 44, 34 


C 


1000-6440 


40 


548 


0, 20, 56, 44 



Table 1: Details of measurements A, B and C. The quoted range of trajectories is the one for 
isrc = and A is the separation in units of trajectories between subsequent measurements for each 
source position t SIC - 



correlation functions obtained using sets B and C were computed with Pn = Pn' =0 and 
the momenta of the initial and/or final pions were induced by twisting one of the pions' 
valence quarks. For each measurement we applied the full twist along one of the spatial 
directions. We changed this direction frequently as the measurements proceeded in order 
to reduce the correlations. In the cases in which both the initial and the final pion carried 
a twist, 6 and 9' were chosen to be anti-parallel. 

Based on a preliminary study of a subset of data set A, we determined the pion mass 
am n to have the central value 0.1907 and this guided us to choose twisting angles so 
as to obtain a suitable range of momentum transfers. (After a detailed analysis by the 
RBC/UKQCD collaboration on their data set of choice, called the FPQ data set in [1], 
the mass was quoted as 0.1915(8).) For such a mass, as mentioned in section 2, the 
minimum value of Q 2 which can be reached without using twisted boundary conditions is 
( a( 2)min ~ 0-051 (Q 2 nin ~ 0.152 GeV 2 ). In order to reach smaller values of Q 2 we introduce 
three twisting angles 2.6832, 2.1285 and 1.6, and in tab. 2 we summarize the corresponding 
kinematics. 

3.3 Three point functions from noise source propagators 

Lattice quark propagators are calculated by inverting the Dirac matrix D upon a matrix 
valued source 77, 

S A ,c(t,x;U) = ^2^2v^ B (t,x;ti,y)i]B,c(y) , (3.6) 

y B 

where A,B,C are spin-colour indices. 

The hadronic form factor calculation is traditionally performed using point source 
propagators [2], for which the Dirac matrix is inverted from a single site with r]B,c(y) = 
8^^ 8b,c- However it has been shown [4,24,25] that the use of stochastic sources allows 
for the calculation of meson propagators at a substantially reduced cost. 

Following [24-26] we use source matrices with random elements from the set Z(2) for 
both real and imaginary components on a single source spin-colour index (0), for all sites 
y on the source timeslice: r]B,o(y, U) £ ^(2) £3 Z(2). With sources of this form, the solution 
S(x, t; ti) requires only a single inversion rather than the 12 required for the point solution. 



data set \p\L \p'\L aE„(p) aE„(p') (aQ) 2 Q 2 (GeV 2 



f wn (q 2 



B 








0.1910(4^ 


0.1910(4) 








1 


B 





1.6 


0.1910(4) 


0.2023(4) 


0.004 


0.013 


0.9804(15) 


B 





2.1285 


0.1910(4) 


0.2106(4) 


0.007 


0.022 


0.9660(24) 


B 





2.6832 


0.1910(4) 


0.2213(4) 


0.012 


0.035 


0.9477(36) 


C 


1.6 


1.6 


0.2023(6) 


0.2023(6) 


0.018 


0.053 


0.9189(75) 


C 


2.1285 


1.6 


0.2106(5) 


0.2023(6) 


0.024 


0.072 


0.8943(88) 


C 


2.1285 


2.1285 


0.2106(5) 


0.2106(5) 


0.031 


0.094 


0.867(10) 


C 


2.6832 


1.6 


0.2213(5) 


0.2023(6) 


0.031 


0.094 


0.864(11) 


C 


2.6832 


2.1285 


0.2213(5) 


0.2106(5) 


0.040 


0.120 


0.838(12) 


c 


2.6832 


2.6832 


0.2213(5) 


0.2213(5) 


0.050 


0.150 


0.802(15) 


A 








0.1912(7) 


0.1912(7) 








1 


A 


2tt 





0.3242(4) 


0.1912(7) 


0.051 


0.152 


0.809(14) 


A 


yj22lt 





0.4167(3) 


0.1912(7) 


0.086 


0.258 


0.711(26) 



Table 2: Table of accessible values of Q 2 = —q 2 for the matrix element (ir(p')\V\ir(p)} together 
with the values of f 7 '' K (q 2 ) ■ For data set B and C we also determined the correlation functions with 
momenta \p\L and \p'\L interchanged. 



A set {rf\j = 1, . . . , N} of these sources has the property that in the limit N — ► oo 



1 N 



(3.7) 



3=0 



such that the pseudoscalar two-point correlator at zero momentum tends to the spatial 
average of the point source solution [24, 25] 



a(t,0)=^^tr{ 7 ^(x,t;t,h 5 (7 5 ^(^*;^)7 5 ) t } 

j=0 x 



(3.8) 



■r-Jl 



Although this explicitly projects to zero momentum at source, twisted boundary conditions 
can be used to induce a non-vanishing meson momentum. 

The properties of equation (3.7) are retained on average when the sources rf reside 
on different configurations such that the stochastic sum can be included in the ensemble 
average. Therefore we require only a single stochastic source per configuration, giving an 
overall factor of 12 cost reduction over the traditional method. 

This technique can be extended simply to three-point correlators using standard se- 
quential source methods 



S'(ti;tf,pf,t,x) = ^75 [V 1 (t,x;t / ,x / )7 5 S'(t / ,x / ;tj)e" 



■ipf-Xf 



t 



75 



(3.9) 



x f 



the solution of which is again non-zero only on a single source spin-colour index, thus 
requiring only one extra inversion. The stochastic cancellation with the other source occurs 
at the source timeslice U as in (3.8). 

3.4 Electromagnetic form factor of a pion with m^ = 330 MeV 

The main results of our computation, the form factor of a pion with m n = 330 MeV for 
a range of low values of Q 2 (obtained from all 3 data sets A, B and C), are presented in 
tab. 2 and plotted in fig. 2. The quoted energies have been determined using the continuum 
dispersion relation E(p) = y/m% + p 2 , where m n is the measured pion mass. We also show 
a zoom of the points at the lowest values of Q 2 . The vertical dashed line represents the 
position of the lowest value of Q 2 which can be reached with periodic boundary conditions 
(Qmin — 0.15 GeV 2 ). From the figure it is satisfying to observe that at Q 2 nin the results 
obtained with partially twisted boundary boundary conditions join smoothly onto the data 
obtained by performing a Fourier sum with momentum of magnitude 2ir/L. 

Our results from datasets B and C are well represented in the range < Q 2 < Q^m 
by the phenomenological pole formula 

&( ^T^V <3 ' 10) 

From the slope of the form factor at q 2 = we obtain the pion's electromagnetic charge 
radius, which is defined by 



,2\ 



d 

'dq~ 2 



(<)-6— rV)L 2=0 . (3.11) 



The best fit, which is shown as the blue curve in fig. 2, gives {r 2 )^30MeV = 6/M 2 ^ = 
0.382(37) (12) (15) fm = 0.382(42) fm , where the first error is statistical and the second is 
due to the uncertainty in the lattice spacing. The third error is to account for our lack of 
a continuum extrapolation (as discussed in sec. 4.2 below). This result corresponds to a 
pole mass of (aM po i c ) 2 = 0.202(20). 

We compare our results to those of the UKQCD/QCDSF collaboration [27] who deter- 
mined the pion form factor for a number of unphysical pion masses m n > 400 MeV using 
periodic boundary conditions. For each pion mass, they fit their data to the pole form in 
(3.10) and hence determine the dependence of the pole mass M po i c on the pion mass. Their 
results are well described by the ansatz, 

M 2 (m 2 ) = c + cim 2 , (3.12) 

for which they determined Co = 0.517(23) GeV and c\ = 0.647(30). Thus, for a pion of 
mass 330 MeV they predict (^) 330 McV = 0.396(15) fm 2 . This result, which we also 

illustrate in fig. 2, is compatible with ours. 

Although the pole formula (3.10) is a good representation of our data for the full 
range Q 2 < Q^mi we ^ n d that the points at the smallest values of Q 2 tend to give a 
slightly smaller central value for the charge radius. We will take as our best estimates 
of (r 2 )33QMcV the value obtained by applying SU(2) ChPT to the points at small Q 2 as 
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Figure 2: Results for the form factor / 7r7r (<7 2 ) for a pion with m, = 330 MeV. The blue clashed 
curve is a pole fit to the data, while the red dashed curve shows the prediction for a 330 MeV pion 
using results from the QCDSF/UKQCD collaboration [27]. The lower plot is a zoom into the very 
low Q 2 region. 



explained in section 4 (we find (r 2 )330MeV = 0.354(31) fm 2 , see tab. 5 ). If we limit the fits 
to the points at small Q 2 , the slope at Q 2 = (and hence the charge radius) is not sensitive 



data set 


maximum Q 2 


linear 


quadratic 


cubic 


pole 


B 

B 
B 


0.013 GeV 2 
0.022 GeV 2 
0.035 GeV 2 


0.354(28)(11) 
0.354(26)(11) 
0.353(25)(11) 


0.353(35)(11) 

0.355(32)(11) 


0.351(41)(11) 


0.361(29)(12) 
0.364(27)(12) 
0.366(27)(12) 


C 


0.150 GeV 2 


0.332(28)(11) 


0.387(44)(13) 


0.406(56)(13) 


0.382(37)(12) 



Table 3: Results for (7^)330 McV obtained by fitting to linear, quadratic or cubic functions of Q 2 
and by using the pole ansatz (3.10). In the first row we use only the single point at the lowest 
value of Q 2 (Q 2 = 0.013 GeV 2 ), in the second we use the two points at the lowest values of Q 2 
(Q 2 = 0.013 GeV 2 and Q 2 = 0.022 GeV 2 ) and in the third row we use the points at the lowest 
three values of Q 2 (Q 2 = 0.013 GeV 2 , Q 2 = 0.022 GeV 2 and Q 2 = 0.035 GeV 2 ). The final row 
corresponds to fits to all 9 points with Q 2 < Q 2 ^^. The two quoted errors are statistical and that 
due to the uncertainty in the lattice spacing. 

to the precise form of the fitting function. To illustrate this we present in tab. 3 the results 
obtained by fitting our results for the form factor at the lowest three values of Q 2 to the 
pole form (3.10) as well as to linear, quadratic and cubic polynomials. In the final row of 
tab. 3 we present the value of (r 2 )330MeV obtained by applying the same fits to all 9 points 
up to Q 2 nin . 

3.5 Comparison of the cost of using point source and Z(2)-wall source propa- 
gators 

In this study we have used two different formulations of the source in the computation of 
the quark propagators. The correlation functions on data set A have been computed from 
point source propagators while the correlation functions on data sets B and C have been 
computed using the noise source technique briefly described in section 3.3. 

In this section we compare the relative computational cost of each approach in order to 
achieve similar statistical errors for standard observables relevant for the phenomenology 
of light mesons. In a very similar recent study [4] such a comparison was carried out for the 
meson spectrum on a 16 3 x 32 lattice also using Nf = 2 + 1 Domain Wall fermions and the 
Iwasaki gauge action. On this smaller volume the inverse lattice spacing was found to be 
a -1 = 1.63(3) GeV and the study was performed using a pion with mass am n ~ 0.44. The 
statistical error on the pseudoscalar and vector meson correlation functions was studied at 
a fixed computational cost, i.e. at a fixed number of inversions of the Dirac matrix. It was 
found that the stochastic (one-end) approach offers a factor of two reduction in the error 
and a definite improvement in plateau quality over the traditional point source technique. 
Preliminary results indicating similar improvements were also reported by ETMC in [3]. 
Here we compare the costs for both approaches on a larger volume and for a much smaller 
pion mass. In particular we perform the comparison for am n , Zy and f 7T7T (q 2 = -Q^m)- 
Table 4 shows the results for each quantity for data sets A and C. In the second column we 
give the number of inversions of the Dirac matrix that were carried out in each case. For 
one measurement 12 inversions are necessary in the case of point source propagators while 
only one inversion is necessary when using the noise source technique. On data set A we 
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data set inversions m^ Zy f W7r {-Qn 



mm / 



A 6444 0.1912(7) 0.7148(9) 0.809(14) 

C 548 0.1910(6) 0.7136(8) 0.802(15) 

Table 4: Comparison of cost and error on quantities relevant for light meson phenomenology. 

have 179 x 12 inversions times three for the number of sources. Our results for data sets A 
and C indicate that the same statistical error for m^, Zy and /"(- Q m j n ) can ^ e achieved 
with only about l/12th of the computational cost when using the noise source technique. 
(This approximate gain of a factor of 12 found for this particular simulation should not be 
confused of course with the 12 inversions performed for each configuration and source for 
data set A.) 

We have also tried to study the error for point-source and noise-source correlators at 
fixed cost, i.e. for a given number of inversions. The cost of the 1176 measurements which 
we carried out with the noise source (data set B) corresponds to 1176/12 = 98 point source 
measurements. While we could carry out reliable fits to the correlators on data set B this 
was not the case for the sub-set of 98 measurements of data set A and no quantitative 
comparison seems possible. This observation shows however that the statistical properties 
of the correlation functions determined with noise-source propagators are better at the 
same computational cost. 

Very light chiral quarks will display near zero modes associated with topological objects 
sampled in the ensemble. Intuitively, we might expect that point source propagators are 
more susceptible to the corresponding fluctuations, particularly if the location of the source 
is in the vicinity of such near zero modes. By contrast we expect such outliers to be averaged 
away when using a volume source like the one considered in this work. Using this picture 
it is not surprising that the gain observed in tab. 4 goes far beyond the one for the error 
of the pion mass observed in [4]. In that work the pion mass was more than double the 
one used here, thus the density of near zero modes was smaller. Furthermore the volume 
was (2/3) 3 of the one used here, allowing for less volume averaging in the case of the noise 
source. 

4. Electromagnetic form factor of a physical pion 

Having determined the electromagnetic form factor of a pion with m w = 330 MeV we 
can estimate what we would expect for that of a physical pion. The natural approach to 
perform this extrapolation is chiral perturbation theory and in the following subsection 
we briefly summarize the predictions of both the SU(2)l x SU(2)r and SU(3)l x SU(3)r 
theories. 

4.1 Chiral perturbation theory for the pion electromagnetic form factor 

The electromagnetic form factor of the pion has been studied extensively in both SU(2)l x 
SU(2)r and SU(3)l x SU(3)r chiral perturbation theory (ChPT). NLO expressions ap- 
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pear in [28, 29] with extensions to NNLO in [30-32] and we now briefly summarise the 
results at NLO. NLO calculations in quenched ChPT and partially-quenched ChPT ap- 
pear in [33]. An NLO calculation with partially twisted boundary conditions in partially- 
quenched ChPT exists in [34]; this is particularly useful to estimate the finite- volume effects. 
In our lattice simulation we use unitary points (each valence quark mass is matched by 
a sea quark mass) and have small finite- volume effects (this will be justified below). We 
therefore use the continuum (unquenched) QCD results to obtain the form factor of the 
physical pion. 

Current conservation ensures that f ww (0) = 1. At NLO only one low energy constant 
(LEC) is relevant for the form factor in both the SU(2) and SU(3) cases. This is denoted by 
Iq(h) for SU(2) and Lg(/i) for SU(3) where the superscript r stands for 'renormalised' and 
we have explicitly indicated the dependence on the renormalization scale /U. The SU(2) [28] 
and SU(3) [29] expressions for the form factor are: 

/su(2),nlo(9 2 ) = 1 + J 2 [-2*6 Q 2 + mml , q 2 ^ 2 )] (4.1) 

/su(3),nlo(9 2 ) = 1 + 4 \ 4L 9 I 2 + 4 «K , q 2 , i, 2 ) + 2H(m 2 K , q 2 ,/, 2 )] (4.2) 



o 



where 



and 



*,( 2 2 2^ m 2 H(q 2 /m 2 ) q 2 m 2 

H{m '« '" > = —3^ 192^ l0? V ( 3) 



H(x) = -J + ^x - ^-^\— - log v ,* ^/^ (4.4) 




with H(x) = —x/6 + 0(x 3 ' 2 ) for small x. For the space-like form factor considered in 
this paper x = q 2 /m 2 is negative and (x — 4)/x > 1 so that the logarithm in (4.4) is real 
as expected. / and /o are the pion decay constants in the SU(2) and SU(3) chiral limits 
respectively (rn u = m^ = with m s at its physical value for SU(2) and m u = m^ = m s = 
for SU(3)). 

The NLO expressions for the charge radius are: 

/ 2\ 12 ^6 1 f, m l A /lr . 

<^)SU(2),NL0 = —£ - ^272 log-f + 1 , (4-5) 



P 8vr 2 / 2 V M 2 



24LS 1 /, ml \ 1 /, m 2 . 



WW** = ~jf - ^ (l° g -f + 1J - J^J! ( log ^ + 1 j ■ (4.6) 

Comparing the expressions for the charge radius gives the relation between the SU(2) and 
SU(3) NLO LEC's [29]: 



in 



IHH) = -2Ll(n) + __^log-f + lj , (4.7) 

where m\ is the kaon mass in the chiral limit for the light quarks. Using the rho-mass 
for the renormalization scale, \i = m p , the second term on the right hand side of this 
relation is very small compared to the expected (power-counting) size of the LECs, so that 
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/g(?TT,p) « —2L r 9 (m p ). A word of caution should be added however. In deriving eq. (4.7) 
from eqs. (4.5) and (4.6) we have set /o = / which is correct at this order. In ref. [1] it was 
found that ///o — 1-23 and so we may expect significant corrections to (4.7). We follow 
the approach of ref. [1] and use SU(2) ChPT to obtain our best results. 

The formulae above are obtained in infinite volume. Jiang and Tiburzi have used 
partially quenched, partially twisted SU(2) chiral perturbation theory to evaluate the finite- 
volume effects in the case where only one of the valence quarks is twisted [34] . This is the 
case for our three points at the lowest values of Q 2 (Q 2 = 0.013 GeV 2 , Q 2 = 0.022 GeV 2 , 
and Q 2 = 0.035 GeV 2 ), which are the points which we use to determine the charge radius 
and the LECs Iq and Lg . From figures 7 and 8 of [34] we see that for the pion mass 
(m n = 330 MeV) and volume, (2.74 fm) 3 , used in our simulation, the finite volume effects 
in (r 2 ) and in 1 — f n1T (q 2 ) are less than 1%. Since the remaining errors quoted for these 
quantities for a pion with m^ = 330 MeV are 7-8%, we feel confident in neglecting the 
finite volume effects in the remainder of this analysis. In order to extend the calculations 
of ref. [34] to the case in which more than one of the valence quarks satisfies twisted 
boundary conditions we would have to perform SU(3) ChPT calculations in the partially 
quenched, partially twisted theory [2]; this is left for a future publication. 

The recent QCDSF/UKQCD two-flavour results for the pion charge radius [27] include 
a larger systematic error, a downwards shift of 6-7%, arising from finite volume effects. The 
estimate of these effects in [27] is obtained using a very different approach to the one we use. 
QCDSF/UKQCD fit form factor data from a range of lattice ensembles, each with a range 
of pion masses, to a pole form, eq. (3.10), with the pole mass given by eq. (3.12). A chirally- 
extrapolated value for the pole mass translates directly to the chirally-extrapolated result 
for the charge radius. Finite volume corrections are modeled by replacing the expression 
for the pole mass with 

M 2 (m 2 ) = c + cim 2 +c 2 e- m - L , (4.8) 

where L is the spatial extent of the lattice. For this second form, additional lattices with 
varying volumes are added to the fit, but the results for the lightest pion, 400 MeV, are 
omitted. A chiral and infinite volume extrapolation now yields a new physical charge 
radius, with the difference quoted as a finite volume systematic error. 

We end this subsection with a discussion of another source of uncertainty which the 
use of chiral perturbation theory can help to estimate. The mass of the (sea) strange quark 
(m s ) in the simulation is different from the physical one (am s = 0.04 in the simulation 
compared to the physical value 0.0343(16) found in ref. [1]). In SU(3) ChPT we use the 
mass of the kaon as found from our simulation and hence obtain the value of the LEC Lg 
without the need for further corrections. The LEC Iq of SU(2) ChPT on the other hand 
depends on the mass of the strange quark and, since this is our preferred approach, we need 
to understand the amount by which Iq could be shifted due to the different value of m s . 
Using eq. (4.7) and the value of the mass from [1] to estimate rhx, we find that the shift 
in ^g(m p ) is about 0.9% and is hence negligible compared even to the 9% statistical error 
(11% total error) found in section 4.2 below (this is also the case if we use eqs. (4.5) and 
(4.6) without setting / = /o, when the relative error grows to 1.3%). For the remainder of 
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the analysis we therefore neglect this uncertainty. 



4.2 Results for the physical pion 

ChPT describes the behaviour of the form factor as a function of both the momentum 
transfer and the quark masses, providing that these are sufficiently small. We fit our data 
at fixed quark masses (i.e. for the pion with mass 330 MeV) as a function of q 2 to the NLO 
formulae for both SU(2) and SU(3) ChPT, eqs. (4.1) and (4.2) respectively. In these fits 
we use the results af = 0.0665(47) and a/o = 0.0541(40) which were determined by the 
RBC/UKQCD collaboration in [1] (in our normalization the decay constant of the physical 
pion is f n ± = 130.7(4) MeV). In this way we obtain the LECs l§ and Lg. Having obtained 
the LECs in this way, we then use the ChPT formulae given above to determine the form 
factor (and hence the charge radius) of a physical pion (m n = 139.57 MeV [35]). 

In ref. [1] it was found that whereas both SU(2) and SU(3) ChPT fit the data for 
the pion masses and decay constants, in the SU(3) case the NLO corrections were very 
large, particularly for the decay constant, casting doubt on the convergence of the chiral 
expansion. For this reason, in ref. [1] the main results were obtained using SU(2) ChPT and 
the above result for the decay constant in the chiral limit, af, includes both the statistical 
and systematic errors. The corresponding result for the decay constant in the SU(3) limit, 
a/o, on the other hand, includes only the statistical error. 

The results of the chiral extrapolation are summarized in tables 5 and 6 for the SU(2) 
and SU(3) cases respectively. In both tables the first column corresponds to the result of 
fitting only to the data point at our lowest value of Q 2 (Q 2 = 0.013 GeV ) to determine 
the single LEC {l§{m p ) or Lg(m p )) and the charge radius. In the second column we use 
the data points at the lowest two values of Q 2 (Q 2 = 0.013 GeV 2 and 0.022 GeV 2 ) and in 
the final column we fit the data for the lowest three values of Q 2 . The results in the three 
columns do not show any dependence on the chosen fit range at these small values of Q 2 . 

Our simulation was performed at a single value of the lattice spacing and we cannot 
extrapolate our results to the continuum limit. However, our action has 0{a 2 ) discretization 
errors and we follow [1] by assigning a systematic uncertainty of 4% to measured quantities, 
representing an estimate of (oAqcd) 2 for our lattice spacing. Thus we assign a 4% error 
from this source to our values for 1 — /^^(o 2 ). This relative error is propagated to our 
results for the LECs and (r 2 ), where it appears as the last error quoted in tables 5 and 6. 

Based on the experience of ref. [1] and because we only know the statistical error for 
a/o, we take for our best estimate the result from the fit to the SU(2)l x SU(2)r expression 
at NLO including the three data points at Q 2 = 0.013, 0.022 and 0.035 GeV 2 , 

l r e (m p ) = -0.0093(10), <r 2 ) 3 30McV = 0.354(31), (r 2 ) x = 0.418(31) . (4.9) 

Comparison of our values for lg(m p ) and Lg(m p ) in tables 5 and 6 with the SU(2)-SU(3) 
conversion formula in (4.7) reveals deviations up to around 50%. By this we mean that the 
LECs obtained directly from the fit differ from the values extracted using the conversion 
formula with the other LEC as input. Large SU(3) NLO corrections were seen in the 
analysis in [1], and indeed the discrepancy can be reduced very significantly by using 
eqs. (4.5) and (4.6) without setting / = /q. 
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Q 



LxIGeV 2 



0.013 



0.022 



0.035 



100 F 6 (m p ) 

(rl)330MeV [fm 2 
<r 2 > x [fm 2 ] 



-0.932(79)(03)(63)(40) 
0.354(28) (12) (00) (14) 
0.418(28)(12)(04)(14) 



-0.933(73)(03)(63)(40) 
0.354(26)(12)(00)(14) 
0.419(26)(12)(04)(14) 



-0.932(71)(03)(63)(40) 
0.354(25) (12) (00) (14) 
0.418(25) (12) (04) (14) 



Table 5: Results from the SU(2) ChPT fits. The errors are statistical, uncertainty in the lattice 
spacing, uncertainty in af and uncertainty from the continuum extrapolation respectively. The 
three columns correspond to using the data at the lowest, the lowest two and the lowest three non- 
zero values of Q 2 respectively, while Q^ax denotes the largest value of Q 2 used in the determination. 



Qi 



'GeV 



0.013 



0.022 



0.035 



100 L r g (m p ) 0.307(26)(03)(49)(13) 0.308(24) (03) (49) (13) 0.308(23)(03)(49)(13) 

(r- 2 >330Mev[fm 2 ] 0.354(28)(12)(00)(14) 0.355(26)(12)(00)(14) 0.355(25)(12)(00)(14) 



(r 



7T/X 



fm 2 



0.460(28)(12)(16)(14) 0.460(26)(12)(16)(14) 0.460(25)(12)(16)(14) 



Table 6: Results from the SU(3) ChPT fits. The errors are statistical, uncertainty in the lattice 
spacing, (statistical) uncertainty in afo and uncertainty from the continuum extrapolation respec- 
tively. The three columns correspond to using the data at the lowest, the lowest two and the lowest 
three non-zero values of Q 2 respectively, while Q^ ax denotes the largest value of Q 2 used in the 
determination. 



collaboration 




technique 


(r£) x [fin 2 ] 


PDG 


[35] 




0.452(11) 


— 


Nam, Kim 


[36] 


instanton vacuum, large N c 


0.455 


. 


QCDSF/UKQCD 


[27] 


N f = 2 Clover 


0.441(19) 


JLQCD 


[37] 


N f = 2 Clover 


0.396(10) 


- 


JLQCD 


[38] 


N f = 2 Overlap 


0.388(15) 


— 


RBC/UKQCD 


this work 


N f = 2 + 1 Domain Wall 


0.418(31) 


— 



Table 7: Previous determinations (excluding quenched lattice results) of the pion's charge radius 
together with the value from the Particle data Group. 

In tab. 7 we compare our result for the charge radius to the one determined from 
experiment and to other recent computations. Note that the previous lattice results were 
obtained with 2 flavours of sea quarks (Nf = 2) and using periodic boundary conditions 
so that the values of Q 2 are much larger than in this paper. 

In fig. 3 we plot our lattice data for the 330 MeV pion and the form factor of a physical 
pion obtained from this data using SU(2) ChPT. The experimental data from ref. [39] is 
also plotted together with the ChPT formula with the PDG world average for the charge 
radius (see also tab. 7). 
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Figure 3: Comparison of experimental results (magenta diamonds) for the form factor f V7r (q 2 ), 
lattice results at m n = 330 MeV (grey triangles and dash-dotted grey line) and the extrapolation 
of the lattice results to the physical point (blue solid line) using NLO SU(2) chiral perturbation 
theory. In addition we also represent the PDG world average for the charge radius using the black 
dashed line. 



5. Summary and conclusions 

In this paper we have successfully used partially twisted boundary conditions to compute 
the electromagnetic form factor of a pion with mass 330 MeV at low values of Q 2 . We use 
our results to compute the LEC l§ of NLO SU(2) and then to determine the physical form 
factor and charge radius, see (1.1). We are able to calculate the form factor for values of 
Q 2 below the minimum value accessible with periodic boundary conditions, see fig. 2. The 
results which we obtain are in good agreement with the experimentally determined form 
factor which gives us further confidence in the use of chiral perturbation theory in the mass 
range below 330 MeV (indeed the value of / which we use in the chiral extrapolation was 
obtained with pion masses up to 420 MeV in ref. [1]). The techniques used in this paper 
can also be applied to other flavour non-singlet form factors of mesons and baryons and 
we strongly advocate the use of partially twisted boundary conditions in order to improve 
significantly the momentum resolution in lattice phenomenology. 

One limitation of the current calculation of the pion's electromagnetic form factor is 
that it was performed at a single value of the lattice spacing, albeit with an action for which 
the discretization errors are of 0(a 2 ) and with good chiral and flavour properties. We are 
currently generating a set of configurations with the same action on a 32 3 lattice with a 
finer lattice spacing and will repeat the present calculation with this ensemble. Although 
the mass and momentum transfers are sufficiently small to expect that NLO SU(2) ChPT 
is a good approximation, it would be nice to be able to check this explicitly. It is not clear 
whether in practice a full NNLO calculation can be performed with sufficient precision (i.e. 
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whether the NNLO LECs will be determined sufficiently accurately) but, as it becomes 
possible to reach lighter quark masses, in the future we will be able to check the stability 
of the results. The finite-volume corrections for our mass and volume are small [34] and 
with our precision can be neglected. 

In our calculation, we confirm the significant reduction in computational cost when 
computing three-point correlation functions using propagators computed from a single 
time-slice stochastic source compared to using point-source propagators. 
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